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ABSTRACT 

The Lyman a emission line is an essential diagnostic tool for probing galaxy formation and evolution. 
Not only is it commonly the strongest observable line from high-rcdshift galaxies but from its shape 
detailed information about its host galaxy can be revealed. However, due to the scattering nature of 
Lya photons increasing their path length in a non-trivial way, if dust is present in the galaxy the line 
may be severely suppressed and its shape altered. In order to interpret observations correctly, it is 
thus of crucial significance to know how much of the emitted light actually escapes the galaxy. 

In the present work, using a combination of high-resolution cosmological hydro-simulations and an 
adaptively refinable Monte Carlo Lya radiative transfer code including an advanced model of dust, 
the escape fractions / esc of Lya radiation from high-redshift {z = 3.6) galaxies are calculated. In 
addition to the average escape fraction, the variation of / esc in different directions and from different 
parts of the galaxies is investigated, as well as the effect on the emergent spectrum. 

Escape fractions from a sample of simulated galaxies of representative physical properties are found 
to decrease for increasing galaxy virial mass M vir , from / csc approaching unity for M v - lr ~ 10 9 M 
to / osc less than 10% for M vir ~ 10 12 M©. In spite of the dust being almost grey, it is found that 
the emergent spectrum is affected non-uniformly, with the escape fraction of photons close to the line 
center being much higher than of those in the wings, thus effectively narrowing the Lya line. 
Subject headings: galaxies: high-rcdshift — radiative transfer — dust — extinction — scattering — 
line: formation — line: profiles 



1. INTRODUCTION 

Many astrophysical and cosmological key questions de- 
pend upon precise measurements of the luminosities of 
distant galaxies; in particular, galactic star formation 
rates (SFRs) and histories, as well as luminosity func- 
tions (LFs), are crucially contingent on the amount of 
assumed luminosities. A fundamental problem in this 
context is naturally the question of how large a fraction 
of the emitted light actually escapes the galaxy. If an 
unknown fraction of the emitted light is absorbed, either 
by gas or by dust, the inferred quantity of interest clearly 
will be subject to large uncertainties or, at best, a lower 
limit. 

The Lya line, albeit notoriously challenging to inter- 
pret, is an extremely powerful probe of the high-rcdshift 
Universe. At z > 2.1, the cosmological redshift of Lya 
makes it the strongest emission line in the optical-NIR 
window, while for z > 4, Lya emitters (LAEs; galaxies 
detected from their emission in Lya, either by narrow- 
band photometry or spectroscopy) becomes easier to 
detect than Lyman-break galaxies (LBGs; galaxies de- 
tected from their dropout in wavelengths blueward of 
the Lyman-break). 

From the shape of the Lya line profile, a variety of 
valuable information can be obtained about the galax- 
ies emitting it. Lya resonant scattering theory as well 
as simulations tell us that, in general, due to the large 
optical depth of the neutral hydrogen for a line center 
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photon, the radiation must diffuse to either the blue or 
the red side of the center, and consequently should es- 
cape the interstellar medium (ISM) of its host galaxy in 
a double-peaked profile. The exact shape of this pro- 
file depends on the physical state of the gas, such as its 
column density, temperature, velocity field, and, in par- 
ticular, dust contents. Similarly, the strength of the line 
depends on the stellar population and the dust contents. 
Notwithstanding the complexity, and indeed sometimes 
intertwinement, of these dependencies, this also makes 
Lya a potentially strong source of information. 

Whereas absorption processes in gas in many cases are 
well-known, the effect of dust on the radiative transfer 
(RT) is still an intensely debated subject. In contrast 
to gas, laboratory experiments with dust are extremely 
complex, partly due to the complications involved in 
replicating the physical environments of the ISM, partly 
due to our limited knowledge vis-a-vis what actually con- 
stitutes cosmic dust. 

In the present-day Universe, most dust is formed 
in the atmospheres of stars on the asymptotic giant 
branch (AGB) of the Hertzsprung-Russel diagram; the 
dying phase of stars less massive than ~8 M© (e.g. 
iHofner fc Andersenll2007HMattsson et al.ll2008l ). In these 
environments the gas is sufficiently cool, yet sufficiently 
dense that molecules may form and stick together to form 
dust grains. 

However, there is observational evidence that dust 
is also abundantly present in the ear ly Universe (e.g. 
IStratta et all 120071 : ICoppin et aT]|2009() . Since the time 
to reach the AGB phase is of the order of 1 Gyr, some- 
thing else must have provided the ISM with dust at these 
epochs. A promising candidate is supernovae (SNe), the 
ejecta of which are believed to exhibit favorable con- 
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ditions for the formation of dust for a short period of 
time, approximate ly 600 days after the explosion (e.g. 
iKottak et aLll2009h . 

For the Milky Way (MW), as well as the Small and 
Large Magellanic Clouds (SMC; LMC), the dust ex- 
tinction curves, i.e. the extinction of lig ht at a given 
wave l ength, are fairly we ll established (e.g.lBianchi et al.l 
fl99l iNandv et all Il98l iPrevot et all 119841 ; iPeil 1199^) 7 
and from the observed color excess E(B — V) one may 
then derive the total extinction. The term "extinction" 
refers to removal of light from the line of sight, be it due 
to absorption or scattering, and may be characterized by 
the number A\ of magnitudes by which the observed ob- 
ject is diminished. For more distant galaxies one is usu- 
ally obliged to assume similar extinction curves. Since 
the stellar population of the SMC is younger than that of 
the LMC, an SMC extinction curve might be expected be 
able to describe better the dust in high-redshift galaxies, 
and has indeed proved to be a g ood fit in GRB host galax- 
ies (e.g. Uakobsson et al.ll2004f ) . Note, however, that the 
prominent feature at 2175 A, characteristic of the LMC 
and MW extincti on, has been detected in a few cases also 
at high redshift (lJunkkarinen et al.l 120041: lEllison et al.l 
120061: ISrianand et al.ll2008t lEhasdottir et al.ll2009D . 

The overall normalization of extinction curves comes 
from the observed property that the extinction is found 
to be very close to pr oportional with the column density 
N-h of hydrogen fe.g. lBohlin et al.lfl978l) . Typically, one 
combines measurements of AThi (and A^h, ) with the ex- 
tinction in the V-band, Ay . In this way, one then knows 
how much light is extinguished when traveling a given 
physical distance in space. 

However, for light that does not travel directly from 
the source to the observer, as is the case for resonantly 
scattered lines like Lya, the situation becomes more rav- 
eled. Not only does the total distance covered by the 
photons increase by a large and a priori unknown factor, 
but the photons received from a given point on the sky 
may also have traveled through physically different en- 
vironments, in turn implying an unknown and possibly 
highly increased probability of being absorbed by dust. 

For this reason, the observed fact that Lya radiation 
nonetheless does escape has long puzzled astronomers. 
The fact that Lya line profiles are often seen to ex- 
hibit a P Cygni-like profile has led to the suggestion that 
high-velocity outflows of gas are needed to enable escape 
(IKunth et al.lll998l : lOstlin et al.l[2008t lAtek et alJl2008h . 
However, at high redshifts many galaxies are still accret- 
ing matter, which should result in an increased blue peak. 
Since this is rarely observed, the shape could be caused 
by other mechanisms, e.g. IGM absorption. 

The angle under which a galaxy is viewed may also 
affect the amount of observed radiation. Ionizing UV 
radiation could create "cones" of low neutral hydro- 
gen density emanating from the star-forming regions 
through which the Lya can escape ([Tenorio-Tagle et al.l 
119991: iMas-Hesse et al.ll2003t) . Even without these ionized 
cones, scattering e ffects alone may cause an anisotr opic 
escape of the Lya (jLaursen fc Sommer-Larsenll2007|) . 

Another commonly repeated scenario is a multi-phase 
medium, where the dust is locked up in cold clouds so 
that the photons primarily travel in an ioni zed, dustless 
medium (Neufeld 1991 1 lHansen~fcO h 200G). Since con- 



tinuum radiation travels through the cloud, it would be 
attenuated more by the dust. This could explain the 
high Lya equivalent widths occasionally observed in in 
LAEs (e.g.lMalhotra fc Rhoadsll2002l ; iRhoads et aTll2003l ; 
iShimasaku et al.ll2006f) ~ 

Previous attempts to determine Lya escape frac- 
tions from high-redshift galaxies have mainly been try- 
ing to match observed Lya luminosities with expected, 
and different methods obt ain quite different results. 
iLe Delliou et~aTl (l2005l[200l found very good agreement 
between galaxies simulated with the galaxy formation 
model GALFORM and observational data at z = 3-6, us- 
ing a constant escape fr action of f e S c = 0. 02 and assum- 
ing no IGM absorption. iDave etaLl (2006) obtained sim- 
ilar results by matching the Lya LF of galax ies from their 
cosmo logical SP H simulation to the dat a of lSantos et al.l 
(|2004h . although iNagamine et all (|2008ft argued that the 
data is based on a small sample and that the simulation 
box size is too small . Matching the simulated Lya LF to 
the ob served one bv lOuchi et alj l|2008[ ). [Nagamine et al.l 
(|2008f ) themselves obtain / esc ~ 0.1, although the pre- 
ferred scenario is not that a certain fraction of the Lya 
radiation escapes, but rather that a certain fraction of 
LAEs are "turned on" at a given time (t he so-called 
"duty cycle scenario"). In a similar way, iDaval et al.l 
(|2009f ) find somewhat higher escape fractions at z ~ 5.7 
and ^6.5 (f eac ~ 0-3), which they use for predicting the 
L F of LAEs at z ~ 7.6. 

iGronwall et al.l (|2007f ) compared inferred Lya and rest- 
frame UV continuu m SFRs of a large sample of LAEs 
from the MUSYC (jGawiser et al.ll2006h survey and ar- 
gue that an escape fractio n of ~l/3 is needed t o explain 
the discrepancy, although iNilsson et all (|2009f) pointed 
out that a missing (1 + z) -factor probably explains the 
difference. Matching Lyq-inferred SFRs to SED model- 
ing of observed LAEs. lGawiser et al.l (|2006f ) found an / esc 
of ~0.8, with a lower limit of 0.2. While SED fitting may 
not be the most accurate way o f estimating SFR s . aim - 
ing to match these observations. iKobavashi et al.l (|2007f ) 
obtain similar result theoretically by incorporating the 
effects of galactic outflows. 

In this paper we aim to scrutinize the effect of dust on 
the Lya R T by applying the M onte Carlo RT code Mo- 
CaLaTA (jLaursen et al.ll2009T) on a number of simulated 
galaxies, extracted from fully cosmological TreeSPH sim- 
ulations. The theoretical aspects of cosmic dust is dis- 
cussed in and the implementation of this into the RT 
code in The new version of the code is then tested 
against physically idealized situations for which analyt- 
ical solutions exist in SjU before it is applied to realistic 
simulated galaxies in Sj5] The results, and the sensitiv- 
ity of these on the values of different input parameters, 
are presented in $6] and Sj7l respectively. Finally, the ob- 
tained results are summarized and discussed in 

2. THEORY OF DUST 

Four quantities characterize what impact the dust 
grains will have on the propagating Lya photons: the 
density] the (wavelength dependent) cross-section of in- 
teraction; the albedo giving the probability that a pho- 
ton incident on a dust grain will be scattered rather than 
absorbed; and finally the phase function defining the di- 
rection into which a non-absorbed photon is scattered. 
These quantities will be discussed below. 
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Dust grains are built up from metals, and thus the dust 
density is expected to scale with gas metallicity in some 
fashion. Metals are created in dying stars, i.e. in AGB 
stars and SNc. For sufficiently dense and cold environ- 
ments, the neutral metals form molecules which eventu- 
ally stick together to form dust. No formal definition of 
the distinction between large molecules and dust grains 
exist, but may be taken to be of the order of ^500 atoms 
or so. 

Depending on the abundances of the individual met- 
als, as well as the physical conditions, a variety of dif- 
ferent types of dust may be produced, with regards to 
both composition and structure, and hence with differ- 
ent scattering properties. Much effort has been put into 
unraveling the nature of cosmic dust, in particular in 
explaining the 2175 A bump. This feature is generally 
attributed to carbonaceous materials, e.g. graphite, di- 
amonds, and/or polycyclic aromatic hydrocarbons, but 
still the precise nature remains unknown. 

In principle, the result of a photon interacting with 
a dust grain may be calculated analytically by solving 
Maxwell's equations, on the basis of the geometry of the 
particle and its optical properties, i.e. the dielectric func- 
tions. This is trivial in the case of simple geometries such 
as spheres and spheroidals. More general shapes and 
composites can be modeled by discretizing the grain into 
a large numb er of dipoles; the so-called discrete dipole ap - 
proximation (jPurcell fc Pennypackerll973l ; lDra ine f988), 
but for the complex and, more importantly, uncertain or 
unknown shape of realistic grains, this is not feasible. 

Had we full knowledge of the relevant properties of 
dust, a distribution of the various species could be cal- 
culated in simulated galaxies, and the radiative transfer 
could then be realized by computing the total optical 
depth of the ISM as a sum of all contributors, and deter- 
mining for each scattering the kind of particle responsi- 
ble for the scattering. However, lacking a sound theory 
of the formation of dust grains, in particular in the high- 
rcdshift Universe, we take a different approach: although 
the exact nature of cosmic dust is not known, the aver- 
age extinction — and hence the cross-sectional area — 
of dust as a function of wavelength is known for many 
different sightlines through the SMC and the LMC (e.g. 
iGordon et alJ 120031 ). Since the metallicity of the Mag- 
ellanic Clouds is fairly well known, the extinction curve 
of the SMC (or LMC) can be scaled to the metallicity 
of the simulated galaxies, thus yielding the extinction in 
the simulations. 

2.1. Cross-section 

Obscrvationally, the extinction Ay in the U-band is 
found to have a surprisingly constant proportionality 
with the column density of hydrogen from sightline to 
sightline within the MW (e .g. iMassa fc Fitzpatrick|[l986l ; 
iFitzpatrick fc Massal 120071 ) . Similar results, but with 
different normalizations, are found for the SMC and 
the LMC (jGordon et al.ll2003h . Accordingly, the cross- 
section crd(A) of dust may be conveniently expressed as an 
effective cross-section per hydrogen atom, thus eliminat- 
ing any assumptions about the size distribution, shape, 
etc., and merely relying on observed extinction curves. 
The optical depth Td of dust when traveling a distance r 
through a region of hydrogen density ??h is then 

Td = n-ara d = A H <Td- (1) 
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Fig. 1. — Extinction cross-section fits to the observed extinction 
curves of the LMC (solid) and the SMC (dashed). The difference 
in amplitude is mainly due to the SMC being less metal-rich than 
the LMC. The vertical, grey-shaded area is the region inside which 
the (rest-frame) Lya line is expected to fall. The two inlets show a 
zoom-in of this region on the extinction curves (top: LMC, bottom: 
SMC), demonstrating the linearity across the Lya line. 

The quantity usually measured is Ax/Na, and the 
cross-section is then 
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We use the fit to the S MC or LMC extinction curves 
proposed bv iPeil (1I992D. wh i ch is an extension of the 
iMattis. Rumpl. fc Nordsieckl (|1977t) -model. The fit is a 
sum of six terms (Drude profiles) representing a back- 
ground, a far-ultraviolet (FUV), and a far-infrared (FIR) 
extinction, as well as the 2175 A, the 9.7 /im, and the 
18 /xm extinction features. Based on newer data from 
IWeingartner fc Draind (|2001h . iGnedin et ail (|2008f ) ad- 
justed the fit and added a seventh term to account for 
the narrow, asymmetric FUV peak in the dust extinc- 
tion. 

Figure Q] shows these fits. The inlets show that the 
extinction curves are very close to being linear in the 
vicinity of the Lya line. In fact, in this region it is an 
excellent approximation to write the cross-section as 
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o- d /10 

0.395 + 1.82 x 10~ 5 (T/10 4 K) 1 / 2 x for the SMC 
0.723 + 4.46 x 10" 5 (T/10 4 K) 1 / 2 x for the LMC 



(3) 



Here, we have used the parametrization x = [y — 
fo)/A^D of the frequency v of the photon, where vq is 
the line center frequency and Ai^d is the Doppler width 
of the line (times ^/2) resulting from the thermal motion 
of the atoms with temperature T. Note that T only en- 
ters Eq. ([3]) to account for the temperature dependency 
of x; Cd itself is independent of T. 

2.2. Number density 

The reason for the variability of extinction with galaxy, 
and the non-variability with sightline, is to a large de- 
gree the different overall metallicities of the galaxies. Al- 
though differences exist within the gala xies, the d iffer- 
ences are larger from galaxy to galaxy (jPeil 119921 ) . In 
most of the calculations we will use an SMC curve, but 
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as shown in $7] the result is not very different if an LMC 
curve is used. 

Because the cross-section is expressed as a cross-section 
per hydrogen atom, the relevant quantity is not dust den- 
sity, but hydrogen density. However, since in general the 
mctallicity at a given location in a simulated galaxy dif- 
fers from that of the Magellanic Clouds, the amplitude of 
the extinction will also differ. Assuming that extinction 
scales with metallicity, a corresponding pseudo number 
density rid of dust at a given location of hydrogen density 
tt-h and metallicity Zj of element i can then be calculated 



J2i Z i,0 ' 



(4) 



where Z^q is the average metallicity of element i in the 
galaxy the extinction curve of which is applied. Obvi- 
ously, rid is not a true dust number density, but merely 
a rescalcd hydrogen number density. 

On average, the SMC metallicities of the different ele- 
ments are deficit relative to Solar values by 0.6 dex (e.g. 
IWeltv et al.lll997T). wh ile the LMC is deficit by 0.3 dex 
(e.g. lWeltv et al.lll999ft . Small metal- to- metal deviations 
from this exist, but scaling Zj to the metallicity of the 
indivi dual metals, using values from iRussell fe Dopital 
(|1992j) . does not change the outcome significantly (cf. $7]). 

The reason that Eq. ([4]) is not expressed as a strict 
equality is that we have so far neglected to differentiate 
between neutral and ionized hydrogen. Dust grains may 
be destroyed in a number of ways, e.g. through collisions 
with other grains, sputtering due to collisions with ions, 
sublimation or evaporati on, or even explo sions due to 
ultraviolet radiation (e.g. lGreenberdll976[ ). These sce- 
narios are all expected to become increasingly important 
for hotter environments. Accordingly, studies of the in- 
terstellar abundances of dust have usually assumed that 
ionized regions contribute negligibly to the dust density, 
and merely concerned themselves with measuring column 
densities of neutral hydrogen, i.e. Hi + H2. Moreover, 
many metallicity measurements are derived from low- 
resolution spectra not capable of resolving and charac- 
terizing various components of the ISM. As discussed in 
g!2.2. 1L dust is also observed in regions that are primarily 
ionized, and since the bulk of the Lya photons is pro- 
duced in the proximity of hot stars with a large intensity 
of ionizing UV radiation, even a little dust associated 
with the ionized gas might affect the results. 

Hence, we assume that the amount of dust scales with 
the total amount of neutral hydrogen plus some fraction 
/i on of the ionized hydrogen, and Eq. ((4| should then be 

V Z 

n-d = (n H i + /ion^Hn) ^ ' - (5) 

Of course, this is not a physical number density of dust 
grains but with this expression, the total optical depth of 
gas and dust as seen by a propagating photon traveling 
a distance r is 

not = r(nmo- x + n d a d ), (6) 
where the neutral hydrogen cross-section a x is given by 



fi 
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Here /12 is the Lya oscillator strength, e and m e is the 
charge and the mass of the electron, and c is the speed 



of light. The wavelength dependence of the cross-section 
enters through the line profile 4>(x) = H (a, x)/^tt, where 
H(a, x) is the Voigt function — the convolution of a 
thermal and a natural broadening of the line — with 
a = Avl /2Avo being the ratio between the natural and 
(twice) the thermal line width (the "damping parame- 
ter"). 

In principle, the summation term in Eq. ([5]) should also 
include a term accounting for the fact that the dust-to- 
metal ratio /dm in a given cell may be different from that 
for which the empirical data exist. In the Milky Way 
and the Magellanic Clouds, /d m — 1 for most metals, 
i.e. roughly 1/2 of the metals is condensed to dust grains. 
The depletion patterns in high-redshift galaxies are not 
well constrained, but no measurements suggest that it 
should be substantially different from the local Universe. 
In fact lPei. Fall, fc Haus er (1998) interpret the depletion 
pat terns of Cr and Zn m easured in damped Lya systems 
by iPettini et all (|1997f ) as giving /d m ~ 1 for z < 3. 
Similarly, fitting de pletion patterns of ei ght elements in 
GRB host galaxies. ISavaglio et all (|2003l ) find /d m — 1. 

2.2.1. Dust in ionized gas 

Ionized gas is found in a number of physically distinct 
locations throughout the Universe. Compact Hn regions, 
or Stromgren spheres, surround young, hot stars, while 
more diffuse Hn is a part of the ISM. Larger H11 "bub- 
bles" are formed around regions of massive star formation 
due not only to ionizing radiation from the stars but also 
to the energy deposited in the ISM from supernova feed- 
back. Outside the galaxies, the IGM is predominantly 
ionized at z < 5-6. Observations show or indicate the 
presence of dust in all of these media. While generally 
lower than in the neutral gas, inferred dust-to-gas mass 
ratios (/dg) in ionized gas span a range from roughly 
equal to the typically assumed MW ISM value of ~0.01, 
to upper limits of ~10~ 4 times lower than this. 

Based on 45- 180 urn (FIR) spectroscopy, 

lAannestad fc Emervl (|2001f ) found the Galactic Hn 
region SI 25 to be strongly depleted of dust, with a 
dust-t o-gas ratio of /d g < 10~ 6 , while ISmith et alJ 
using MIR imaging and spectroscopy, inferred 
a dust-to-gas ratio of the Galactic Hn region RCW 
38 of 10" 5 t o 10~ 4 . O n the other hand, using FIR 
spectroscopy IChini et"aTI (|1986f ) found 12 Hn regions 
to be dust-depleted by "only" a factor of 10 relative 
to the MW ISM (i.e. ,f dg ~ 10~ 3 ), while from FIR 
photometry, lHarper fc Low! (|l97lh found the median 
dust-to-ionized-gas ratio of seven Hn regions to be close 
to 0.01. 

For the more diffuse Hn gas that comprises part of the 
ISM, most obtained extinction curves in a sense already 
include the contribution of Hn to nd, although its quan- 
tity is not revealed when measuring Hi column densities. 
Hence, any value of /i on for the diffuse ISM larger than 
would account twice for the ionized gas. 

The dominant destruction mechanism of dust is prob- 
ably shock waves, a ssociated with, e.g., high - velocity 
clouds and SN winds (jDraine fc SalDetedll979m . How- 
ever, since SNe are thought to be the prime creator of 
dust at high redshifts, the Hn bubbles in the vicinity 
of massive star-forming regions cannot be entirely void 
of dust, and observational evidence of dust related to 
SN remnants (SNR) and starburst regions does indeed 
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exist. Using MIR imaging. iBouchet et ail (|2006fl deter- 
mined the dust-to-gas ratio of SN 1987A to be ~5 x 10~ 3 . 
Somewhat lower re sults are found in Kes 75 (~10~ 3 
from FIR and X-rav. iMorton et al.ll2007t ) and in Kepler 's 
SN (~1(T 3 from IR and bremsstrahlung,|CmrtinJ|200j). 
On larger scales, the hostile environments imposed by 
the SNe and ionizing radiation will reduce the dust 
density in starbur s t regi ons. Fitting continuum SEDs, 
iContini fc COTrtiml (|2?)0l found that 10~ 4 < f dg < 1CT 2 
in various starburst regions in a sample of seven luminous 
infrared galaxies. However, such regions are not ionized 
to the same level as compact Hii regions and SNRs, and 
as argued in the case of the diffuse Hn, the scaling of dust 
with Hi to some extend already accounts for the Hn. 

Various feedback processes are also responsible for ex- 
pelling a non-vanishing amount of metals and dust into 
the IGM, although inferred dust-to-g as ratios ten d to be 
small: from IR-to- X-ray luminosites, iGiard et al.l (|2008f ) 
inferre d a dust-to-gas r atio o f a few to 5 times 1CT 4 , 
as did IChelouche et alJ (|2007l ) by comparing photomet- 
ric and spectroscopic properties of quasars behind SDSS 
clusters. Higher (d ust-to-Hi ~ 0.05 in the M81 Group, 
IXilouris et al.ll200a l — possibly expelled from the star- 
burst galaxy M82 — and lower (/dg ~ 10~ 6 in the 
Coma cluster and even less in five other Abell clusters, 
IStickel et al.ll200"2l ) values are also found. Additionally, 
sputtering by the hot halo gas may tend to destroy pri- 
marily small grains, leading to a flattening of the ex- 
tinction curve in the UV; at the Lya wavelength, this 
may reduce the aver age cross-section by a factor of 4-5 
(jAguirre et al.ll2001l ). 

In summary, the factor f- lon is a practical way of mod- 
eling the destruction of dust in physically "hostile" envi- 
ronments. For simplicity, in the RT code we will not dis- 
tinguish between Hn in various regions but merely settle 
on an average dust-to-gas ratio of ionized gas of ~ 10~ 4 ; 
that is we set /i on = 0.01. In Sj7l we investigate the effect 
of other values of /; on and find that using 0.01, the result- 
ing escape fractions lie approximately midway between 
those found when using /i on = (corresponding to the 
complete destruction of dust in regions where hydrogen 
is ionized) and /i on = 1 (corresponding to no destruction 
of dust at all). Moreover, these extreme values does not 
seem to change / esc by more than ~ 25%. 

2.3. Albedo 

When a photon interacts with a dust grain, it may be 
either absorbed or scattered. The efficiency with which 
the dust grain scatters radiation is dependent on the 
composition (material, shape, etc.) of the dust and on 
the wavelength of the incident photon. If the photon 
is not scattered (i.e. emitted with the same wavelength 
as the incident photon), it is absorbed. In this case it 
is converted into heat and re-emitted at a later time as 
infrared radiation. Expressing the total cross-section as 
a sum of a scattering cross-section a s and an absorbing 
cross-section er a , such that ad = a s + er a , the albedo A of 
the dust is defined as 



A = ^. 

0~d 



(8) 



wavelength, A lies approximately between 0.3 and 0.4 for 
various size distributions fitted to the LMC and SMC, as- 
sumin g that the dust is made mainly of graph ite and sili- 
cates |Peilll992t IWeingartner fe Drainej|2q01l). W e adopt 
an albedo of A = 0.32 (from lLi fc Draindl2001h . but in- 
vestigate the impact of using other values in SjJl 

2.4. Phase function 

If a photon is not absorbed, it is scattered. The prob- 
ability distribution of deflection angles 6 from its orig- 
inal path is given by the phase function. For reasons 
of symmetry, the scattering must be symmetric in the 
azimuthal angle cf> (unless the grains are collectively ori- 
ented in some preferred direction due to, e.g., magnetic 
field lines), but in general this is not the case in 8. In fact, 
dust is often observed to be co nsiderably forward scatter- 
ing [e.g. in re flection nebulae (iBurgh et al.ll2002l ). diffuse 
galactic light (|S~chiminovich et al.ll2001h . and interstellar 
clouds dWitt fc Oliverilll990ft ]. This asymmetric scatter - 
ing may be described by the iHenvev-GreensteTnl (|1941[ ) 
phase function 



^hg(m) = = - 



1-9 2 



The albedo of dust has been investigated observa tion- 
ally from reflection nebula e fe.g.lCalzetti et al.lll995[ ) and 
diffuse galactic light (e.g. iLillie et al.lll976f h~At the Lya 



2 (l+ 5 2_ 2^)3/2' (9) 

where fi = cos 9 and g = (fj,) is the asymmetry param- 
eter. For g = 0, Eq. © reduces to isotropic scattering, 
while g = 1 (—1) implies complete forward (backward) 
scattering, g is a function of wavel ength, but for A close 
to that of Lya, iLi fc Draini (|2001[ ) found that g = 0.73. 
Again, other values are investigated in [J7] 

3. IMPLEMENTATION IN MONTE CARLO CODE 

The numerical code, MoCaLaTA, used to conduct 
the radiative t ransfe r of Lya, is described in detail in 
lLaursen et al.l (|2009f ). To understand how the effect of 
dust is implemented, a brief summary of the concepts of 
the code is given below. 

The non-dusty version of the code assumes an arbitrary 
distribution of neutral hydrogen density tihi, Lya emis- 
sivity L, gas temperature T and velocity field Vbuik- The 
physical parameters, typically resulting from a cosmo- 
logical simulation, are assigned in a cell-based structure, 
where any cell may be refined, i.e. split up in eight sub- 
cells, recursively to an arbitrary level of refinement, thus 
allowing for investigation of a detailed structure. 

The basic principles of the code are as follows: pho- 
tons arc emitted from a random point in space, weighted 
according to the emissivity at that particular point, i.e. a 
given photon has a probability Li/L to t of being emitted 
from the i'th cell, where L to t is the total Lya luminosity 
of the system. It is emitted in a random direction, with 
a frequency close to the Lya line center. 

The photon travels through the gas, traversing an opti- 
cal depth t with a probability e~ r , before it is scattered 
on a hydrogen atom. Re-emitted in a new, random direc- 
tion, weighted according to an appropriate phase func- 
tion, it continues its journey until it escapes the compu- 
tational box. At each scattering, the velocity (thermal 
+ bulk) of the atom will Doppler shift the frequency of 
the photon, so that it diffuses not only in real but also in 
frequency space. This procedure is then repeated for a 
number n p h of photons sufficient to provide good statistic 
for the quantities of interest (several 1000s, or millions, 
depending on the desired output). 
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For each photon and for each scattering, the probabil- 
ity that the photon escape in the directions of six virtual 
observers located at a distance of the luminosity distance 
of the galaxy along the three Cartesian axes is computed 
and added to an associated array (a "CCD") of three di- 
mensions; two spatial and one frequential, so that a full 
spectrum is obtained for each pixel in the image. 

3.1. Metallicity normalization 

In addition to the parameters riHij L, T and Vbuik; 
every cell has an ionized hydrogen density rtHn and a 
metallicity Zi of the different elements associated with 
it, from which a dust density (per hydrogen atom) n d is 
calculated according to Eq. ([5]). 

With a randomly drawn optical depth to be covered by 
a photon, the distance r traveled can then be calculated 



(10) 



Having reached this distance, a random number 1Z G 
[0, 1] (a "univariate") determines whether the photon hits 
a hydrogen atom or a dust grain by comparing it to the 
ratio g = n<i<r<i/ (nm<rx + «dO"d) = Td/(r x + r d ); if TZ < g, 
the interaction is caused by dust. In this case a second 
univariate is compared to the albedo of the dust grain, 
dictating whether the the photon is absorbed, thus termi- 
nating the journey of this particular photon, or scattered, 
in which case it is re-emitted in a random direction given 
by Eq. ©. 

Various acceleration schemes were implemented in the 
non-dusty Lya RT; the extension of these to the dusty 
version is discussed in App. 

4. TESTING THE CODE 

iNeufeldl (|1990f ) provided an analytical expression for 
the escape fraction of photons emitted from inside a 
"slab" (i.e. finite in the z-direction and infinite in the xy- 
direction) of an absorbing medium. The solution, which 
is valid for very high optical depths (aro > 10 3 , where To 
is the optical depth of neutral hydrogen from the center 
to the surface of the slab) and in the limit (aro) 1 / 3 3> T a , 
where r a is the absorption optical depth of dust, is 



cosh CV(oto) 1/3 



(11) 



where C' = a/3/C 71 " 5 ^ 12 ^ w hh ( ~ 0.525 a fitting parame- 
ter. Figure [2] shows the result of a series of such simula- 
tions, compared to the analytical solution. 

Various AMR grid configurations were tested. Of 
course, the physical parameters such as njji, n d , and T 
of a cell does not depend on the refinement level, but the 
acceleration schemes discussed in App. lAldo. 

5. SIMULATIONS 

To investigate the fraction of emitted Lya pho- 
tons that escapes galaxies at high redshift, the code 
was applied to a number of different simulated high- 
resolution galaxies, emerging from a fully cosmologi- 
cal N-body/hydrodynamical Tree SPH simulation. This 
code was described in det a il in | Sommer-Larsen et al.l 
(l200l and IS ommer-Larsenl (|2006f ) . and summarized in 
lLaursen et al.l ( 2009f ). Here, we will only briefly review 
the basic characteristics of the code: 
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Fig. 2. — Escape fractions / csc of photons emitted from the 
center of a semi-infinite slab of gas damping parameter a, hydrogen 
optical depth to, and dust absor bin g optical depth r a , compared 
to the analytical solution in Eq. mil l. 

The numerical simulations arc first carried out at low 
resolution, but in a large volume of space. Subsequently, 
interesting galaxy-forming regions are re-simulated at 
high resolution. In addition to hydrogen and helium, 
the code follows the chemical evolution of C, N, O, Mg, 
Si, S, Ca, and Fe. 

The Lya emissio n is produced by three differen t 
processes (see also lLaursen fe Sommer-Larser] |2007|) . 
viz. from recombinations in photoionized regions around 
massive stars (responsible for ^90% of the total Lya lu- 
minosity), gravitational cooling of infalling gas (~10%), 
and a metagalactic UV background (UVB) photoionizing 
the external parts of the galaxy (~1%). 

The UVB field is a ssumed to be that given by 
lHaardt fe Madaul (|1996l ) , where the gas is treated as op- 
tically thin to the UV radiation until the mean free path 
of a UV photon at the Lyman limit becomes less than 0.1 
kpc, at which point the gas is treated as optically thick 
and the UV field is "switched off" . 

A far more realistic UV RT scheme was implemented 
in the code bv lRazoumov fe Sommer-Larsenl(|2006l . l2007f) 
which alters the ionization state and temperature field 
of the gas somewha t. However, as was concluded in 
lLaursen et "all (|2009Q . the effect of the improved UV RT 
only results in minor changes for the non-dusty Lya RT. 
Moreover, as shown in |j7]it does not alter / osc drastically, 
and has thus not been applied in the present paper. 

Nine individual galaxies are extracted from the cosmo- 
logical simulation at redshift z = 3.6 — at which time 
the Universe was 1.8 Gyr old — to be used for the Lya 
RT. These galaxies are representative of typical galaxies 
in the sense that they span three orders of magnitudes in 
mass, the most massive eventually evolving into a disk 
galaxy with circular speed V c ~ 300 km s" 1 at z = 0. 
The numerical and physical properties of these galaxies 
are listed in Tab. Q] and Tab. [2] respectively. 

Since the Lya RT code is cell-based, the relevant phys- 
ical properties of the particles are interpolated from the 
50 nearest neighboring particles onto a grid of base reso- 
lution 128 3 cells. Dense cells are subdivided in eight cells, 
which are further refined until no cell contains more than 
ten particles. Achieving the same resolution without this 
adaptively refined mesh would require from ^16 000 3 to 
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Characteristic quantities of the simulations 








Galaxy 


S33sc 


K15 


S29 


K33 


S115 


S87 


S108 


S115sc 


S108sc 


JV p ,tot 


1.2xl0 6 


2.2x10 s 


1.1x10 s 


1.2x10 s 


1.3x10 s 


1.4x10 s 


1.3x10 s 


1.3x10 s 


1.3x10 s 




5.5xl0 5 


1.0x10 s 


5.1x10 s 


5.5x10 s 


6.4x10 s 


7.0x10 s 


6.3x10 s 


6.4 x10 s 


6.3x10 s 


mgp H ,m s tar 


5.4xl0 5 


9.3xl0 4 


9.3xl0 4 


9.3xl0 4 


l.lxlO 4 


1.2xl0 4 


1.2xl0 4 


2.6xl0 :! 


1.5xl0 3 


m DM 


3.0xl0 6 


5.2x10 s 


5.2x10 s 


5.2x10 s 


6.6xl0 4 


6.5xl0 4 


6.5xl0 4 


1.4xl0 4 


8.1xl0 3 


eSPH,£star 


344 


191 


191 


191 


96 


96 


96 


58 


48 


£ DM 


612 


340 


340 


340 


170 


170 


170 


102 


85 




18 


10 


10 


10 


5 


5 


5 


3 


2.5 



Note. — Total number of particles (JV p tot)i number of SPH particles only (-/Vgpjj), masses (m), gravity 
softening lengths (e), and minimum smoothing lengths (i m i n ) of dark matter (DM), gas (SPH), and star particles 
used in the simulations. Masses are measured in h^ 1 MQ, distances in h _1 pc. 



TABLE 2 

Physical properties of the simulated galaxies 



Galaxy 


S33sc 


K15 


S29 


K33 


S115 


S87 


S108 


S115sc 


S108sc 


SFR/Mq yr" 1 


70 


16 


13 


13 


0.5 


0.46 


1.62 


3.7xl0" 3 


1.7xl0" s 


M„/Mq 


3.4xl0 10 


1.3xl0 10 


6.0 x 10 9 


6.5xl0 9 


2.5x10 s 


1.8x10 s 


4.9x10 s 


2.0xl0 7 


5.9x10 s 


M vil /M Q 


7.6xlO n 


2.8xlO n 


1.7 x 10 11 


1.3X10 11 


2.5xl0 10 


2.1xl0 10 


2.6xl0 9 


4.9x10 s 


3.3x10 s 


r vir /kpc 


63 


45 


39 


35 


20 


19 


10 


12 


5 


[O/H] 


-0.08 


-0.30 


-0.28 


-0.40 


-1.22 


-1.28 


-0.51 


-1.54 


-1.64 


Vc(z = 0)/km s" 1 


300 


245 


205 


180 


125 


132 


131 


50 


35 


-^Lyc/erg s" 1 


1.6xl0 44 


4.5xl0 4S 


2.9 xlO 43 


2.5xl0 4S 


1.3xl0 42 


1.1 xlO 42 


2.6 xlO 42 


4.9 xlO 40 


3.2xlO ss 


i«/,uv/erg s _1 Hz -1 


5.0xl0 29 


6.7xl0 28 


9.3xl0 28 


5.5X10 28 


3.6xl0 27 


3.3xl0 27 


1.2xl0 28 


2.6xl0 25 


1.2xl0 2s 



Note. — Star formation rates (SFRs), stellar masses (M t ), virial masses (M v j r ), virial radii (r v j r ), metallicities ([O/H]), circular 
velocites (V c ), Lya luminosities (LLya)i antl UV luminosities (Z/^uv) f° r the simulated galaxies. All quoted values correspond to 



if no particularly strong winds are present. 

6. RESULTS 

In the discussion of the results, we will first focus on 
one particular galaxy, observed from one direction. We 
arbitrarily chose K15, which at z = becomes an M31- 
sized disk galaxy, and we chose the observer placed to- 
ward negative values on the z-axis (the z_-direction), 
the direction in which most radiation escapes. In £16.41 
we proceed to discuss the extend to which the results are 
representative for other galaxies and directions. 

6.1. Where are the photons absorbed? 

Figure 2] shows the surface brightness (SB) maps of the 
galaxy K15. The left panel shows how the galaxy would 
look if the gas were dust-free, while the right panel dis- 
plays the more realistic case of a dusty medium. Com- 
paring the two images, it is seen that the regions that 
are affected the most by dust are the most luminous re- 
gions. This is even more evident in Fig. [5l where the 
azimuthally averaged profiles of the SB maps are shown. 

The reason for this is two-fold: the most luminous re- 
gions are the regions where the stars reside. Because dust 
is produced by stars, this is also where most of the dust 
is. Since the stars, in turn, are born in regions of high 
hydrogen density, the RT is here associated with numer- 
ous scatterings, severely increasing the path length of the 
photons, and consequently increasing further the proba- 
bility of being absorbed. 



a redshift of z = 3.6, except V c which is given for z = 0. 

~65 000 3 cells; infeasible with present-day's computer fa- 
cilities. 

This level of resolution is rather crucial to the escape 
of the Lya photons. As is evident from Fig. [31 the bulk 
of the photons is produced in regions of very high column 
densites, exceeding AThi — 10 23 cm -2 . Even for metallici- 
ties of only 1/100 Zq, Eq. (jTTJ) implies an escape fraction 
of -10~ 5 (for T ~ 10 4 K). However, Eq. ([Til assumes a 
homogeneous medium. As was ar gued bvlNeufeldl (119911) 
and investigated numerically by lHansen fc Oh! (|2006f ). 
in a multi-phase medium, the Lya photons may escape 
more easily. In the (academic) case of all the dust resid- 
ing in cool, dense clouds of neutral hydrogen which, in 
turn, are dispersed in a hot, empty medium, the Lya es- 
cape fraction may approach unity. The reason is that the 
photons will scatter off of the surface before penetrating 
substantially into the clouds, thus effectively confining 
their journey to the dustless intercloud medium. 

Although this scenario is obviously very idealized, the 
presence of an inhomogeneous medium undoubtedly re- 
duces the effective optical depth, and has indeed been 
invo ked to explain unusually large Lya equivalent widths 
fe.g. lFinkerstein et~ al. 2008). Nevertheless, it is not clear 
to which degree the escape fraction of Lya will actually 
be affected. Other scenarios have been proposed to ex- 
plain the apparent paradoxal escape of Lya, e.g. galactic 
superwinds, as discussed in the introduction. Here we 
present the first calculations of / csc , based on fully cos- 
mological, numerical galaxy formation models, demon- 
strating that generally of the order 5%-30% of the Lya 
radiation escapes the galaxies at redshifts z ~ 3-4, even 
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log(N H1 /cm 2 ) log(lL/erg s 1 cm" 2 ) 

13.4 15.4 17.5 19.5 21.6 23.6 -8.8 -6.9 -5.1 -3.2 -1.4 0.5 




x/kpc x/kpc 



Fig. 3. — Neutral hydrogen column density (-/Vjjj) map (left) and integrated source Lya emissivity (EL) map (right) of the simulated 
galaxy K15. The vast majority of the p hotons are seen to be emitted in the very dense environments. According to the analytical solution 
for the Lyo escape fraction provided by Ncufcld ( 1990, Eq. {TTJ in this paper), virtually all photons should be absorbed by dust, but taking 
into account the dumpiness of the ISM allows for much higher escape fractions. 
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Fig. 4. — Surface brightness maps of the galaxy K15, as viewed from the negative z-direction, without dust (left) and with dust (right). 
Including dust in the radiative transfer affects primarily the most luminous regions. 



The SB in the less luminous regions also decreases. 
Although dust also resides here, having been expelled by 
the feedback of starbursts, only little absorption actually 
takes place here. The decrease in luminosity is mainly 
due to the photons being absorbed in the high density 
regions that would otherwise have escaped the inner re- 
gions and subsequently been scattered by the circumam- 
bient neutral gas into the direction of the observer. This 
is also discernible from Fig. [6] that displays an image of 
the absorbed photons. 

Figure [7] shows the source Lya emissivity of the pho- 
tons that are eventually absorbed, compared to that of 
the photons that eventually escape. Here it becomes ev- 
ident that virtually all the absorbed photons are emitted 
from the central parts, while photons that escape are 
emitted from everywhere. In particular, the radiation 
produced through gravitational cooling escapes more or 
less freely. 



6.2. Effect on the emergent spectrum 

Due to the high opacity of the gas for a Lya photon at 
line center, photons generally diffuse in frequency to ei- 
ther the red or the blue side in order to facilitate escape. 
Consequently, the spectrum of the radiation escaping a 
dustless medium is characterized by a double-peaked pro- 
file. The broadening of th e wings is dictated by the prod- 
uct otq (jHarringtonlll973f ). i.e. low temperatures and, in 
particular, high densities force the photons to diffuse far 
from line center. Since such conditions are typical of 
the regions where the bulk of the photons is absorbed, 
the emergent spectrum of a dusty medium is severely 
narrowed, although the double-peaked feature persists. 
Figure [8] displays the spatially integrated spectra of the 
dustless and the dusty version of K15. 

This interesting result shows that even though the dust 
is effectively grey (the small wavelength dependence of 
Eq. ([3]) does not produce substantially different results 
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Fig. 5. — Surface brightness (SB) profile of the galaxy K15, 
again with dust (solid) and without dust (dotted). Left ordinate 
axis gives the SB as measured at the source while right ordinate 
gives the values measured by an observer at a distance given by 
the luminosity distance of the galaxy. The decrease in SB in the 
less luminous regions is noticable. However, this decrease is for the 
most part not due to photons being absorbed in the hot and ten- 
uous circumgalactic medium but rather reflects a lack of photons 
that in the case of no dust would have escaped the luminous regions 
and subsequently scattered on neutral hydrogen in the direction of 
the observer. 
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Fig. 6. — Image of the locations of absorption of the Lya radi- 
ation that does not escape K15. Effectively, this image shows the 
column density of dust. 

from using a completely wavelength independent cross- 
section), the emergent spectrum is affected in a highly 
"non-grey" fashion: whereas the escape fraction of the 
inner part of the spectrum is of the order 50%, it rapidly 
drops when moving away from the line center. 

6.3. Escape fraction 

In the z_-direction of K15, the fraction of Lya photons 
escaping is 0.14. As mentioned earlier, the z_-direction is 
the direction into which most radiation escapes — with- 
out dust, this direction is ~ 3 times as luminous as the 
^--direction, which is where the least radiation is emit- 
ted. Including dust, because the brightest regions are 
affected the most, the ratio between the luminosity in 



the least and the most luminous directions is somewhat 
reduced, although z_ is still more than twice as bright 
as X-. 

The sky-averaged escape fraction for K15 is 0.16. Re- 
call that in these simulations SMC dust has been applied. 
However, as discussed in S}7] using LMC dust does not al- 
ter the results significantly. 

6.4. General results 

The results found in the previous sections turn out to 
be quite illustrative of the general outcome of Lya RT in 
a dusty medium: photons are absorbed primarily in the 
dense, luminous regions, leading to a reduced luminosity 
in these parts of the galaxies and effectively smoothing 
out prominent features Furthermore, the wings of the 
spectrum experience a strong cut-off. 

6.4.1. Anisotropic escape of Lya 

In general, the radiation does not escape isotropically; 
the ratio of luminosities observed from different direc- 
tions ranges from ~1.5 to ~4. Without dust, these ratios 
are somewhat higher, up to a factor of ~6. Although 
"bright directions" are affected more by the dust than 
less bright directions, the variation in / esc as a function 
of direction is not large, and not very different from the 
sky-averaged / esc - 

6.4.2. Correlation of f esc with galactic mass 

Figure [9] shows the escape fractions of the galaxies as 
a function of the virial masses of the galaxies. Despite a 
large scatter, and although the sample is too small to say 
anything definite, the figure indicates that / esc decreases 
with increasing mass of the host galaxy. In addition to 
the nine high-resoluti on galaxies [of w hich two a ppear in 
two "v ersions" ; with a lSalpeten (|l955f ) and with a lKroupal 
(|1998f ) initial mass function (IMF)], 17 galaxies from a 
low-resolution (Kroupa) simulation are shown. These 
galaxies were chosen according to the criterion that the 
number of star particles is >1000 to ensure an accept- 
able resolution. Similar re sults are found for the escape 
of ion izing UV radiation (jRazoumov fc Sommer-Larsenl 
2009}. The reason is probably a combination of two 
mechanisms: small galaxies have a lower metallicity and 
hence less dust than large galaxies. Furthermore, due to 
their smaller gravitational potential the stellar feedback 
will "puff up" small galaxies relatively more and make 
them less ordered, thus reducing the column density of 
both dust and gas. 

6.4.3. Narrowing of the spectrum 

Figure [TU] shows the emergent spectra of the studied 
galaxies, arranged according to virial mass. In general, 
the same trend as for K15 is seen for all galaxies: the dust 
primarily affects the wings of the profiles. As discussed 
in £16.2) the reason is that the wings of the Lya profile 
are comprised by photons originating in the very dense 
regions of the galaxy having to diffuse far from the line 
center in order to escape, and since this is also where the 
most of the dust is residing, these photons have a larger 
probability of being absorbed. For the two least mas- 
sive galaxies, S108sc and S115sc, the ISM is neither very 
dense nor very metal-rich, meaning that photons escape 
rather easily. Consequently, the line is not particularly 
broadened and most of the photons escape the galaxies. 
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Fig. 7. — Surface brightness maps of the source Lya emissivity, of the photons that are eventually absorbed (left), and the photons 
that eventually escape (right). These images clearly show that absorbed and escaping photons do not in general probe the same physical 
domains. 
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Fig. 8. — Observed spectral distribution of the radiation escap- 
ing the galaxy K15 in the negative 2-direction, with (solid) and 
without (dotted) dust. The vertical, dashed line marks the Lya 
line center. Although the dust is close to being grey, the spectrum 
is not affected in the same way at all wavelengths; the inner parts 
are only diminished by a factor of ~2, while the wings are severely 
reduced. The reason is that the photons in the wings are the ones 
produced in the dense and dusty regions of the galaxy, so these 
photons have a higher probability of being absorbed. 



6.4.4. Extended surface brightness profile 

The fact that more photons are absorbed in the 
bright regions tend to "smooth out" the SB pro- 
files of the galaxies. Young galaxies have often been 
found to be more extended on the sky when ob- 
serve d in Lya, compared to continuum band observa- 
tion (iMoller fc Warrenl Il998t iFvnbo et all [20Q]1 12003 
iNilsson et al.ll2007ft . lLaursen fc Sommer-Larsenl (|2007f ) 
found that, even without dust, resonant scattering it- 
self may cause an extended Lya SB profile. Including 
dust merely adds to this phenomenon, since steep parts 
of the profile are flattened. Figure [Tl] shows this effect 
for three arbitrarily chosen galaxies and directions. Ex- 
tended emission is seen also on larger scales (tens of kpc), 
although this is more likely to be caused by cold accretion 



Fig. 9. — Escape fraction / csc as a function of galactic virial 
mass M v j r . Although the plot exhibits a large scatter, there is a 
clear trend of / eS c decreasing with increasing M v j r . 

onto dark matter halos ([Nilsson et al. 2006). 

6.5. Temporal fluctuations 

Many factors play a role in determining the exact value 
of / esc . Although the metallicity, and hence the state of 
matureness of the galaxy, as well of the size of the galaxy 
seem to be the most significant property regulating / 0S c, 
less systematic factors like the specific configuration of 
gas and stars are also likely to have a large influence. The 
scatter in Fig. [5] is probably due to this effect. To get an 
idea of the fluctuations of / esc with time, RT calculations 
was run for snapshots of K15 from 100 Myr before (z = 
3.8) to 100 Myr after (z — 3.5) the one already explored 
at z = 3.6. In this relatively short time, neither Z nor 
Mvir should not evolve much, and thus any change should 
be due to stochastic scatter. 

Figure [T^] shows the variation of / osc during this time 
interval, demonstrating that the temporal dispersion is 
of the order 10% on a 10 Myr scale. This suggests that 
the scatter seen in f CS c(M v i T ) (Fig. \§§ reflects galaxy- 
to-galaxy variations rather than / csc of the individual 
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Fig. 10. — Emergent spectra of the studied galaxies, without dust (dotted) and with dust (solid), ordered after decreasing virial mass. In 
order to use the same ordinate axis for a given row, some intensities have been multiplied a factor indicated under the name of the galaxy. 
Generally, lines that are broadened by resonant scattering tend to be severely narrowed when including dust. 
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Fig. 11. — Lycr surface brightness (SB) profiles of three of the simulated galaxies as observed from three arbitrary directions, with (solid) 
and without (dashed) dust. While the intrinsic SB profiles are shown in grey, the black lines show the profiles convolved with a Gaussian 
kernel corresponding to a seeing of 0'.'5. In general, the inclusion of dust tends to "smooth out" the profiles, effectively resulting in a more 
extended SB profile. 
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Fig. 12. — Escape fraction / esc from the galaxy K15 as a function 
of time t over a period of 200 Myr, normalized to / csc at t — 1.8 
Gyr. The dispersion of / osc over time is quite small, indicating the 
the dispersion in Fig.[9]is due to galaxy-to-galaxy variations rather 
than the escape fractions of the individual galaxies fluctuating. 
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Fig. 13. — Relative escape fraction / esc from the galaxy K15 as 
a "function" of model, i.e. simulation where all parameters but one 
are equal to th at o f the benchmark model used for the simulations 
in Sec. [5] See ^7. II for an explanation of abcissa labels. Except for 
the factor /i on , the chosen benchmark model appears quite robust 
to varying other parameters. 

galaxies fluctuating strongly with time. 

7. PARAMETER STUDY 

The adopted model of dust clearly involves a multitude 
of assumptions, some more reasonable than others. To 
inspect the dependency of the outcome on the values of 
the parameters a series of simulations of K15 was run, 
varying the below discussed values. The resulting escape 
fractions compared to that of the "benchmark" model, 
used for the simulations in Sec. OH are shown in Fig. [T3l 

7.1. Varied parameters 

The albedo A of the dust grains. The chosen value of 
0.32 is bracketed by the values 0.5 and 0, i.e. somewhat 
more reflective and completely black, respectively. As 
expected, the higher the albedo, the higher the escape 
fraction, but note that even completely black dust re- 
duces / GSC by less than 10%. This is because the bulk of 
the photons is absorbed in the very dense environments, 



where scattering off of one grain in many cases just post- 
pones the absorption to another grain. 

The scattering asymmetry parameter g. The three 
cases g = 0, 1, —1 are tested, corresponding to isotropic 
scattering, total forward scattering, and total backscat- 
tering. The difference from the benchmark model value 
of 0.73 is virtually non-existent; the fact that most of the 
scatterings take place in the dense environments makes 
the transfer be dominated b y scattering on hydrogen. 

The dust cross-section ovi. iFitzpatrick fc Massal (|2007f ) 
showed that the variance of the extinction curves (in the 
MW, normalized to Ay) is approximately 20% at the 
Lya wavelength. Decreasing (increasing) the dust cross- 
section by this quantity increases (decreases) the escape 
fraction as expected, but not by more than ~±5% per- 
cent. Also, a constant cross-section with o& = <7dU=o 
was tested, but with no notable effect. 

Likewise, a variance of the reference metallicity Z 
is present from sightlinc to sightlinc in the Magellanic 
Clouds, probably at least 0.1 dex. Using a smaller 
(larger) reference Z makes the metallicity in the sim- 
ulations comparatively larger (smaller), with a larger 
(smaller) dust density as a result and hence a smaller 
(larger) escape fraction. Since the escaping photons rep- 
resent different sightlines in the galaxies, it is fair to use 
the average Z (for the individual metals) in Eq. (J5]), but 
to investigate the sensitivity on the reference metallicity, 
simulations with Z^mc increased (reduced) by 0.1 dex 
was run, resulting in a —5% (+5%) change in / osc . 

Letting scale with the metallicity of only C and Si 
instead of the total metallicity does not alter / esc much 
either. This model is relevant since C and Si probably 
are the main constituents of dust. 

Although the metallicity of the Magellanic Clouds is 
smaller than that of the MW, the relative abundances be- 
tween the various elements are more or less equal. Small 
deviations do exist, however, but letting n c \ scale with 
the metallicity of the individual metals does not change 
the results much (point labeled "Zindiv" m Fig. [T3)l . 

Dust type. Using an LMC extinction curve instead of 
SMC results in a slightly (few %) larger escape fraction, 
since the quantity cx n^a^ oc <j&/Zq is roughly 10% 
lower for the LMC than for the SMC. 

The minimum temperature T m i„ of the simulations. 
The cosmological simulation includes cooling of the gas 
to ^10 4 K. Since the temperature affects the RT of Lya, 
we artificially lowered the temperature of the cells with 
T ~ 10 4 K to 10 3 , 10 2 (approximately the temperature of 
the cold neutral medium), and 10 K (approximately the 
temperature of a molecular cloud) , to see if not including 
sufficient cooling could affect the results. However, as is 
seen in Fig. 1131 the difference is insignificant. 

The ionizing UV RT scheme. As mentioned already, 
implementing a more realistic UV RT scheme does not 
alter the outcome of the non-dusty Lya RT significantly. 
When including dust, as seen from Fig. 1 131 the improved 
RT results in a slightly increased / osc , although less than 
10%. The reason is that this RT is more efficient than 
the "old" RT scheme at ionizing the neutral gas in the 
immediate vicinity of the stars and, accordingly, at low- 
ering the dust density. However, in these regions the gas 
density is so high that ionization in most cases is followed 
by instantaneous recombination, and hence the physical 
state of the gas in the case of the improved RT is not 



Lya RT with dust: escape fractions from simulated high-z galaxies 



13 



2.0 



S 1.0 



- 


I I I I 1 I 
• 

• 

• 




- 




• 

• 

• 

1 // 1 1 1 1 1 


• 


• 


10~ 3 10" 2 10" 1 1 



fion 

Fig. 14. — Relative variation in escape fraction / eac from the 
galaxy K15, as a function of /i on , the fraction of ionized hydrogen 
that contributes to the density of dust [Eq. JEJ] . Note the discon- 
tinuity on the abcissa axis. Even a little dust in the ionized region 
can affect / esc quite a lot. 

altered significantly. 

The initial mass function. Since a Salpeter IMF is 
more top-heavy than a Kroupa IMF, i.e. produces rela- 
tively more massive stars per stellar mass, it also yields 
a higher metallicity and hence a higher absorption by 
dust. However, the increased feedback from the massive 
stars serves as to counteract star formation, and these 
two effects more or less balances each other. The ratio of 
the Kroupa-to-Salpeter feedback energy is 0.617, while 
for the yield the ratio is 0.575 (for oxygen). The result, 
as is seen from Fig. Q2J is only a slightly smaller escape 
fraction. 

The fraction of ionized hydrogen f lon contributing to 
the dust density. Insufficient knowledge about the dust 
contents of ionized gas is by far the greatest source of 
uncertainty in / esc , as is seen in the figure. This effect is 
further investigated in Fig. [TU where the relative escape 
fraction from K15 for different values of /i on is shown. 

From the figure, it is seen that even a little amount 
of dust associated with the ionized hydrogen can affect 
the escape fraction quite a lot. The reason is that most 
of the scatterings and the absorption take place in the 
dense region where also the star formation is high. In 
these regions supernova feedback Shockwaves will recur- 
rently sweep through the ISM, heating and ionizing the 
medium without significantly lowering the density. For 
/ion = 1, the calculated dust density of these regions is 
not affected, while for / ion = 0, this effect renders the 
gas virtually dustless, resulting in highly porous medium 
with multiple possibilities for the photons of scattering 
their way out of the dense regions. 

The resulting escape fraction of the benchmark model 
lies approximately midway between the two extrema and 
if nothing else, / csc can be regarded as having an uncer- 
tainty given by the result of these extrema, i.e. ^20%. 
Nevertheless, we believe that /i on = 0.01 is a realistic 
value, cf. the discussion in § 12.2.11 

7.2. Resolution dependency 

The galaxies studied have been extracted from a 
medium resolution cosmological simulation and then re- 
simulated at 8 times higher resolution. To check if 



the resolution is sufficient to trust the results, / csc - 
calculations were also performed at the medium resolu- 
tion. The resulting SB profile and spectrum for K15 is 
seen in Fig. [15] The resulting escape fraction is a tiny bit 
lower than for the hi-res simulation, but less than 5%. 

To investigate the significance of the AMR structure, 
we also carried out simulations progressively dcsolving 
the structure level by level. In K15, K33, and S115, the 
maximum level C of refinement is 7, where C = corre- 
sponds to the unrefined base grid of 128 3 cells. Eight cells 
of C = I are desolved to £ = t — 1 by taking the average 
of the physical parameters. Since temperature reflects 
the internal energy of a body of gas, and since the com- 
bined velocity is given by momentum, T and v DU ik arc 
weighted by mass. 

Figure [16] shows the resulting escape fractions of these 
three galaxies. Desolving the first few levels does not al- 
ter / esc notably, indicating that the galaxies are suffiently 
resolved. However, eventually we see the importance 
of the AMR structure: with insufficient resolution, the 
dumpiness of the central, luminous ISM is lost, "smooth- 
ing out" the low-density paths that facilitate escape, and 
consequently f csc decreases. When the resolution be- 
comes even worse, the central regions are averaged with 
the surrounding low-density gas, so that most of the pho- 
tons are being emitted from medium dense cells, resulting 
in a small probability of scattering on neutral hydrogen, 
and hence a small probability of being absorbed by dust. 

8. SUMMARY AND DISCUSSION 

We have undertaken numerical radiative transfer cal- 
culations for Lya in high-redshift galaxies including the 
effects of dust, with special emphasis on the fraction of 
photons escaping the galaxies as well as the impact on 
the emergent spectra. The abundance of, and the inter- 
action probability with, dust is modeled by scaling known 
extinction curves to the (location-specific) metallicity of 
the galaxies, and destruction processes are modeled by 
reducing the dust density in regions where hydrogen is 
ionized. 

The main results are listed below: 

• The escape fraction of Lya seems to decrease with 
increasing virial mass of the host galaxy. This is 
indicated by Fig. [9] which show that / esc is close 
to unity for galaxies of M v ; r ~ 10 9 -10 10 M®, while 
it falls off to a few % for M vil - 10 n -10 12 M Q . 
This effect could be caused partly by the fact that 
smaller galaxies experience less total star forma- 
tion, so that the total amount of metals, and hence 
dust, is smaller than for massive galaxies. However, 
since the same tren d is seen for the escape of ion- 
izing UV radiation (]Razoumov fc Sommer-Larsenl 
|2009( ). for which dust is less important, the chief 
factor is likely to be feedback energy of the smaller 
galaxies being able to "puff up" and disorder its 
host to a larger degree than for massive galaxies, 
due to their gravitational potential being smaller. 

• Although the cross-section of the dust is nearly in- 
dependent of the wavelength of light over the Lya 
line, the spectrum is affected in a highly wavelength 
dependent fashion: close to the line center, the es- 
cape fraction is of the order 50%, while it quickly 
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Fig. 15. — Comparison of the surface brightness profile (left) and emergent spectrum (right) of the radiation escaping the galaxy K15 
when performed on galaxies simulated at low (solid) and high (dashed) resolution. Grey lines show the true profile and black lines show 
the profile convolved with a seeing of C/.'5. While the SB in the central regions appear to agree nicely, the fact that we are comparing two 
different simulations make luminous regions in the outskirt appear somewhat shifted. 

infalling gas. Since the dust originates from stars, 
by far most of the dust is found where the wings are 
produced, meaning that this part of the spectrum 
is affected the most by the dust. 

• For the same reason, the surface brightness profiles 
of galaxies, when observed in Lya, is "smoothed 
out" . Without dust, the central parts of the galax- 
ies are often very luminous compared to the outer 
parts, thus resulting in a central bump in the SB 
profile. With dust, this bump is reduced, giving 
rise to an even more "flat" SB profile, which can 
then be interpreted as an extended SB profile when 
comparing to continuum bands. 

lLaursen fc Sommer-Larsenl (|2007f ) found that even 
without dust, scattering effects alone can explain 
extended Lya profiles, although the very central 
parts may exhibit a steeper profile. Including dust 
in the calculation tends to remove these central 
bumps and thus flatten the profile. 

The obtained results (/ csc , SB, spectra) seem to be 
quite insensitive to the assumed values of various param- 
eters characterizing the dust, such as dust albedo, scat- 
tering asymmetry, dust cross-section, extinction curve, 
etc. Of the studied input parameters, the only actual 
uncertainty comes from insufficient knowledge about the 
dust contents of ionized gas, but this probably at most 
introduces an error of ^20%. This robustness against 
input parameters is convenient in the sense that we can 
rely on the results, but is also a nuisance since at least 
from Lya observations, we should not expect to be able 
to learn much about the physical properties of the dust 
itself. 

As discussed in the introduction, previous attempts to 
determine Lya escape fractions have been quite diver- 
gent, ranging from a few percent to close to a hundred 
percent. In this work we have shown that / esc may in- 
deed vary rather much from galaxy to galaxy, and we 
propose a tentative evidence for a negative correlation 
with galactic mass. Obviously, many factors play a role 
in regulating / csc ; in particular the age of a given galaxy 
will be significant, since the dust accumulates over time. 



Fig. 16. — Relative escape fractions / eac from the three galaxies 
K15, K33, and S115 as a function of the maximum level £ max of 
AMR refinement. £ m ax = corresponds to having only the 128 3 
base grid, while £ m ax < corresponds to desolving the base grid. 
For increasingly lower resolution, / osc drops due to the low-density 
paths being smeared out with high-density regions. Eventually, 
however, when the resolution is so course that the central star- 
and gas-rich regions are mixed with the surrounding low-density 
region, / eS c increases rapidly. 

approaches zero in the wings. Consequently, the 
line is severely narrowed, although its width may 
still reach several 100s of km s _1 . 

The reason is that different parts of the spectrum 
originates in physically distinct environments of 
its host galaxy. The bulk of the Lya photons is 
produced in the central, dense regions, where hot 
stars ionize the surrounding hydrogen which sub- 
sequently recombines, eventually resulting in ^2 
Lya photons for every 3 ionizing photon. Since the 
optical depth of neutral hydrogen is so immense, 
in order to escape the photons have to diffuse in 
frequency to either the blue or the red side of the 
line center, where the cross-section falls off rapidly. 
Hence, the wings of the spectrum emanates from 
the star-forming regions, whereas the central parts 
to a high degree stems from the regions less pop- 
ulated by stars, and from gravitational cooling of 
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Various authors have invoked different scenarios to ex- 
plain their inferred escape fractions, e.g. galactic out- 
flows, ionized paths, multi-phase medium, and viewing 
angle. We find no evidence that any single of these sce- 
narios should be dominating entirely the magnitude of 
/ osc . Rather, a mixture of gas kinematics, ISM dumpi- 
ness and ionization state, as well as viewing angle will 
influence the total observed Lya luminosity, and hence 
the deduced escape fraction. The investigated galaxies 
exhibit both outflows and infall, but not at exceptionally 
high velocites, and artificially setting all velocites to zero 
still allows plenty of radiation to escape (this was done 
for K15; the shape of the spectrum is altered, but / osc 
still is around 0.16). Artificially erasing the dumpiness of 
the ISM decreases / e sc as expected but also in this case 
much radiation still escapes. The luminosity observed 
from different directions may vary by a factor of ~1.5 
to ~4. Note, however, that this factor could increase if 
the effect of AGN was implemented in the cosmological 
simulations. 

We have also investigated the variation in f csc over a 
time span of 200 Myr, and found that / esc exhibits only 



minor fluctuations. 

We have focused on the escape of Lya photons from 
galaxies, and did not take into account any absorp- 
tion/scattering in the IGM. At z = 3.6, the Universe 
is expected to be highly ionized, and thus most of the 
radiation should be able to propagate freely once the 
host galaxy has been escaped. Preliminary RT calcula- 
tions using MoCaLaTA on cosmological volumes indi- 
cate that, on average, approximately 80% of all of the ra- 
diation blucward of the Lya line center should be trans- 
mitted, while the Universe is effectively transparent to 
radiation rcdward of the line center (Laursen et al., in 
prep.). 
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APPENDIX 
ACCELERATION SCHEMES 
Core-skipping scheme 

For a photon with frequency close to the line center, the opacity of the gas is so high that the photon does not 
diffuse significantly in space. In order to accelerate to execution time of the code, these scatterings can be skipped by 
drawing the velocity of the scattering atom not from a pure Gaussian but from a Gaussian where the inner part is 
truncated, such that higher velocities and hence large frequency shifts are favored. This acceleration scheme is invoked 
whenever \x\ <■ x cr jt , where a; cr it varies according to the physical conditions in the cell (for a detailed explanation, see 
lLaursen et afll2009l ). 

With a dusty medium, we must investigate the possibility that the photon would have been destroyed, had we 
not used this acceleration scheme, i.e. the probability -P a bs(^crit) of absorption for a photon initially in the core, 
before escaping the frequency interval [— cc C rit, £c C rit]- Ultimately, we will determine this probability numerically, but 
to interpret the result, we will first investigate the scenario analytically. In the following calculation, factors of order 
unity will be omitted. 

The probability per interaction that a photon with frequency x be absorbed is 



Pabs(£) 



Td 



1 + 0(a;)r o /r a 



(Al) 



since r a ~ Td and t x ~ <P(x)tq. Here, T aj d,2:,o corresponds to the optical depth of this particular part of the journey 
The number dN(x) of scatterings taking place when the frequency of the photon is close to x is 



dN(x) = N tot (j>(x)dx, 



(A2) 



where Ntot is the total number of scatterings before the photon exits [— x C rit, Xcmt], i.e. the total number of scatterings 
skipped. Here we have assumed complete redistribution of the frequency, i.e. there is no correlation between the 
frequency of the photon before an d after the scatteri ng event. This is a fair approximation in the core (IUnnolll952bl ; 
Uefferies fc Whitel ll960'). However. lOsterbrockl (|1962f ) showed that once the photon is in the wing, it has a tendency 
to stay there, only slowly drifting toward the line center with a mean shift per scattering (Ax) = — l/|a?|. 

For the purpose of the current calculation, the Voigt profile is approximated by a Gaussian in the core and a power 
law in the wing, such that 

~ x for x < x cw 



<j>(x) 



for x > 



(A3) 



where x cw marks the value of x at the transition from core to wing (Eq. (15) in lLaursen etafll2009t ). 
At each scattering, the probability of escaping the region confined by a; cr it is 



for x crit < x c 
for a; cr it > x c 



Pcsc(Zcrit) =2 / 4>{x)dx 

( erfc x cr it 

~< _a_ 



(A4) 
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Fig. 17. — Probability of absorption P a b s before escaping the region of the line confined by the value x cr it, for a series of different values 
of aro/ra (labeled at the corresponding lines). 

where crfc is the complimentary error function. 
Using Eq. (|A2[) , the total probability of being absorbed can be calculated as 



Pabs (^crit ) 



Pahs (x)dN(x). = N u 



p & h s {x)(j)(x)dx(x) 



The t otal num ber o f scat terin gs before escape if the photon is not absorbed is -/V to t 
Eqs. l[Xl]) and (|A"4]l . Eq. (|A"5)) then evaluates to 



1/Pe 



For x r 



Pabs (^crit ) 



(7.r 



erfexcrit Jo e x + T Q /r a 



(A5) 



from 



(A6) 



The exponential integral and the factor l/erfcx cr it is of the same order, but since the factor To/r a is of the order 
10 s Z/Zq, Eq. l[A"6]) will usually be negligible. 
In the case of x cr it > x cw , 

f x " il dx 

Pabs(* crit ) ~ / 7T7 vT /-• ( A? ) 

J a/0(x) + flT /T a 

This integral can be evaluated separately for the intervals [0, x cw [ and [x cw , a; C rit]- For the first, the result is usually 
negligible, as in the case with Eq. (|A6j) . The second integral yields 



-fabs (^crit ) 



1 r, _l -^crit 

-(tan 

V t 



where 



. tan" 1 ■), (A8) 

t= (aro/ra) 1 / 2 . (A9) 

For x cr ;t > x cw , the assumption of complete redistribution becomes very inaccurate, as the photon spends consid- 
erably more time in the wings, with a larger probability of being destroyed. However, Eqs. (|A8[) and (|A9[) reveals a 
signature of the behavior of Pabs(a^crit), viz. that it has a "tan~ 1 -ish" shape, and that it scales not with the individual 
parameters a, tq, and r a , but with their interrelationship as given by the parameter t. To know exactly the probability 
of absorption, a series of Monte Carlo simulations were carried out for a grid of different temperatures, gas densities, 
and dust densities. The results, which are stored as a look-up table, are shown in Fig. |T7j Indeed, the same fit applies 
approximately to different T, riHu and giving equal values of t. 

Whenever the acceleration scheme is applied, a bilinear interpolation over log(t) and x cr i t determines the appropriate 
value of P a bs- A univariate is then drawn and compared to P a bs, thus determining if the photon is absorbed or allowed 
to continue its journey. Note that under most physically realistic conditions, only low values of P a bs are actually met. 
However, when invoking the acceleration scheme many times, the probability of absorption may become significant. 

Semi- analytical scheme 

In very dense regions, defined in the simulations as cells with aro > 2000, th e code can be furthe r accelerated by 
calculating analytically the characteristics of a photon escaping such a region. lLaursen et al.l (|2009f ) found that the 
spectrum of photons escaping a sufficiently dense cube is characterized by the analytical solution for the spectrum from 
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Fig. 18. — Escape fractions (black crosses) of photons emitted from the center of a cube of damping parameter a, line center optical 
depth ro, an d dust absorption optical depth T a , measured from the center to the face, compared to the analytical solution (solid black) 
given by Eq. lETot . For comparison, the equivalent results for the slab (grey; the same as in Fig. [2]l are also displayed. 



an equivalent "slab" of gas found bv iHarringtonl (|1973f ) . with the exception that the independent variable otq (where 
ro is measured from the ce nter to the face of th e cube) be replaced by 77070, where n = 0.71 is a fitting parameter (for 
a detailed explanation, see lLaursen et al1l2009f ). 

With the inclusion of dust, we must calculate the possibility of the photon being absorbed in such a cube. From the 
above, we might expect that replacing aro by ijaro and r a by r\r & in the slab-relevant equation for the escape fraction 
[Eq. (fTT|) ] yields the relevant solution. In fact, even better fits can be achieved by also replacing the square root by 
an exponantiation to the power of 0.55. That is, every time the semi-analytical acceleration scheme is invoked, a 
univariate is drawn and compared to the quantity 



/esc 



coshfC [77 4/3 (ar ) 1 /3(i_A)r d 



0.55 



(A10) 



determining whether or not the photon should continue its journey. 

Figure [18] shows the calculated escape fractions from a number of cubes of different physical properties. The tests 
performed in Sj4] and many of the results of S}5] were performed both with and without this acceleration scheme, all 
agreeing to a few percent within statistical errors. 

Luminosity boosting scheme 

Since the vast majority of the photons are emitted within a relatively small volume of the total computational 
domain, many photons are needed to reach good statistics in the outer regions. In order to reach convergence faster, 
the probability of emitting photons from low- luminosity cells can be artificially boosted by some factor l/w > 1, later 
corrected for by letting the emitted photon only contribute with a weight w to statistics (spectra, SB profiles, escape 
fractions) . 

This factor is calculated for the i'th cell as 



W; 



10 log(i i /L max )/6 



(AH) 



where Li is the original luminosity of the cell, L max is the luminosity of the most luminous cell (not to be confused 
with £ ma x), and b is a "boost buffer" factor that determines the magnitude of the boost; for 6 = 1, all cells will have 
an equal probability of emitting a photon while for b — > oo the probability approaches to original probability. 

The optimal value of b depends on the quantity and physical region of interest. Since photons are absorbed primarily 
in the central regions, / osc -calculation will usually converge fastest using b — > oo, and since after all most photons are 
recieved from this region as well, the same counts for the spatially integrated spectrum. If one wishes to investigate 
the SB or the spectrum of the outer regions, b should not simply be set equal to unity, however, since a significant 
fraction of the photons received from here are photons originating in the central parts and later being scattered in the 
direction of the observer. In this case, faster convergence can be reached with b ~ 1.5 to a few tens. 
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